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A fully symmetry unrestricted Time-Dependent Density Functional Theory extended to include 
pairing correlations is used to calculate properties of the isovector giant dipole resonances of the 
deformed open-shell nuclei ^^^Yb (axially deformed), ^**Os (triaxially deformed), and ^^®U (axially 
deformed), and to demonstrate good agreement with experimental data on nuclear photo-absorption 
cross-sections for two different Skyrme force parametrizations of the energy density functional: SkP 
and SLy4. 
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The isovector giant dipole resonance (GDR) is perhaps 
the simplest example of a nuclear collective motion of all 
the protons against all the neutrons. Since its obser- 
vation in the photo-absorption cross section pj, it has 
been intensively studied as it combines several challeng- 
ing aspects of the physics of the atomic nucleus [H 
Even though GDR is practically harmonic in character, 
it is not an adiabatic collective mode and various damp- 
ing mechanisms of the collective energy are at work 0]. 
In the models of Migdal [5], Goldhaber-Teller [6j,and 
Steinwedel-Jensen [7] GDR is described as the relative 
motion of two fluids, either compressible or incompress- 
ible, with neutrons and protons vibrating around a com- 
mon center of mass, and the mass dependence of the ex- 
citation energy reads A~^/^ and A"^/^ respectively [8]. 
A good estimation of the GDR vibrational frequency is 
fiw « 80 MeV for spherical nuclei. GDR is in- 

terpreted simply as the equivalent of the zero-sound in 
a nuclear system and the size of the nucleus sets a con- 
straint on the largest wavelength. In the case of deformed 
nuclei, the GDR peak is split, with various frequencies 
revealing different principal axes of the nuclear shape. 
Since the GDR state is not an eigenstate of the nuclear 
Hamiltonian, it is characterized by a spreading width 
that cannot be extracted from a simple hydrodynami- 
cal approach. The spreading width in principle can be 
calculated within a microscopic model of the atomic nu- 
cleus. The total width of the GDR is mainly due to a 
couple of mechanisms: the coupling of the GDR to com- 
plex nuclear configurations T^, and the coupling to the 
continuum, leading to the escape of neutrons and protons 

. These two widths contribute to the total width of the 
GDR, r = + r^, and their relative contributions vary 
depending on the mass number A and the N/Z ratio. 
The escape width is typically more important for light 
nuclei. The physical mechanisms related to may be 
quite complicated and involve coupling to low energy sur- 
face vibrations, Landau damping and collisional damping 
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The theoretical approach to GDR is typically within 
the linear response of the nuclear system to an external 
probe |9j . In contemporary approaches, the description 
of the atomic nucleus is provided by a Density Functional 
Theory (DFT) 10^ and the GDR is described within the 
small amplitude limit of the time-dependent version of 
DFT. This approach leads to the well known (Quasipar- 
ticle) Random Phase Approximation (Q)RPA that has 
been used with a variety of nuclear density functionals, 
though only very recently for deformed systems, see Refs. 
[T1-21J for several representative calculations. When ap- 
plied to open-shell and, in particular, deformed nuclei 
QRPA requires diagonalizations of matrices of extremely 
large sizes and often severe truncations of the quasipar- 
ticle basis set. Truncations lead to a number of unde- 
sired features, in particular spurious states, instead of 
having a zero excitation energy, get admixed to transi- 
tions of physical interest. In a deformed open-shell nu- 
cleus, states corresponding to the excitation of the cen- 
ter of mass motion, rotation of the system, and gauge 
transformations arising from the breaking of the proton 
and neutron number conservations, all have zero-energy 
excitation energy (sometimes referred to as Goldstone 
modes). When (Q)RPA is derived, one usually allows 
for small violations of the Pauli principle, which have 
only a negligible effect when one is considering the most 
collective transitions. (Q)RPA formally is an approach 
in which the fermionic nuclear Hamiltonian is bosonized 
and the number of bosonic excitations allowed to exist is 
significantly much larger than the number of true phys- 
ical excitations. Significant difficulties arise when one 
tries to go beyond the (Q)RPA and consider coupling to 
more complex configurations, and describing the spread- 
ing width r^^. In spite of the efforts of many generations 
of theorists and many approaches suggested so far (nu- 
clear field theory [22] , boson expansion methods O |23] , 
and diagrammatic methods Pl), a truly satisfactory mi- 
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croscopic treatment is still lacking and theorists resort to 
either phenomenological models, such as optical poten- 
tials (which can easily lead to overestimates of unless 
vertex corrections are accounted for as well or sim- 
ple prescriptions such as adding a popular, but arbitrary, 
energy smoothing of the (Q)RPA strength. 

Here we present an approach to GDR in open-shell nu- 
clei based on an extension of the DFT to superfluid sys- 
tems, superfluid local density approximation (SLDA) and 
its further extension to time-dependent phenomena (TD- 
SLDA), which was developed over the years and applied 
to a number physical systems and phenomena, see Refs. 
[24-35 . The time-dependent Hartree-Fock-Bogoliubov 
(TDHFB) approximation has been used previously as an 
alternative to QRPA in spherical nuclei j36) and with 
a number of further approximations in deformed nuclei 
[371 |3S]. TDSLDA appears formally as a time-dependent 
self-consistent local mean-field approximation. No spuri- 
ous modes are admixed to physical transitions, since the 
energy density functional respects all the required sym- 
metries (translational and rotational symmetry, Galilean 
invariance |3S], parity, isospin symmetry, gauge sym- 
metry) and no further approximations are introduced, 
apart from discretization errors which are under control. 
Isospin symmetry is broken only by the proton/neutron 
mass difference and by the Coulomb energy, which we 
treat in the Hartree approximation. We use an approxi- 
mation to the normal part of the nuclear energy density 
functional provided by various Skyrme force parametriza- 
tions [401 - 142] and all the terms are taken into account 
numerically exactly. The pairing part of the energy den- 
sity functional is treated as described in Ref. [Mj and 
[25j . when a single coupling constant is used for both 
protons and neutrons, and even and odd particle num- 
bers as well, unlike most phenomenological approaches to 
nuclear pairing which break isospin invariance in the pair- 
ing channel. In such approaches, see Refs. [T2] and [TF 
for examples of QRPA calculations, the proton pairing 
coupling constant is larger in magnitude than the neu- 
tron pairing coupling constant, \Vp\ > \Vn\. This leads 
to a flagrant violation of isospin invariance in the energy 
density functional, the magnitude of which is clearly not 
due to charge-symmetry-breaking forces. Such effective 
pairing couplings would be consistent with a new kind 
of pairing energy, that does not violate charge symmetry 
breaking [43^. The Coulomb interaction will lead to just 
an opposite relation [44], namely to \ Vp\ < \Vn\. Gahlean 
invariance [39| requires the presence of currents in the en- 
ergy density functional, and although their contribution 
is vanishing in the ground states of even-even nuclei, it 
is important for the excited states [45l [4^ . 

Within TDSLDA one describes accurately the interac- 
tion in both particle-hole and particle-particle channels, 
the treatment is fully self-consistent, and all symmetries 
of the Hamiltonian are properly accounted for. We place 
the nuclear system on a 3D spatial lattice, needed spatial 



derivatives are determined using Fast Fourier Transform, 
and we at first determine the ground state properties 
within SLDA and subsequently subject the system to an 
external time-dependent one-body potential [531 133 113 ■ 
The Coulomb interaction is treated with particular care, 
so as not to introduce the potential due to image charges, 
which appear naturally in a spatial lattice formulation 
with periodic boundary conditions. Although we study 
the small amplitude limit, the equations we solve here 
are the same as for a motion of arbitrary amplitude. The 
emerging equations are formally equivalent to the TD- 
HFB approximation with local potentials, or to the time 
dependent Bogoliubov-de Gennes (TDBdG) equations: 
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where we have suppressed the spatial r and time coordi- 
nate and k is the label of each quasiparticle wave function 
[uka{'r,t),Vkaii^,t)]. where a =t,i- The single-particle 
Hamiltonian h^^-' (r, t) is a partial differential operator 
(thus local) and A{r,t) is a pairing field, all defined 
through the normal, anomalous, spin, and isospin den- 
sities and currents. The interaction with various applied 
external fields (spin, position and/or time-dependent) is 
described by including the corresponding potentials in 
the single-particle Hamiltonian /icro-' (r, t). This approach 
represents a flexible tool to describe in general large am- 
plitude nuclear motion as it contains the coupling to the 
continuum and between single-particle and collective de- 
grees of freedom, since the meanfield is explicitly time- 
dependent. However, the later type of coupling between 
single-particle and collective degrees of freedom will ac- 
count only for a part of the diagrams discussed in Ref. 
[4]. If the external field will act along one of the sym- 
metry axes for example, only collective oscillations along 
that axis will be excited, but not in the perpendicular 
direction. Although the TDSLDA is designed to provide 
the values of one-body densities, a simple modification 
of the variational principle extends the approach to two 
body observables, for example, particle number fluctu- 
ation |48 - 52j . Another straightforward extension of the 
formalism leads to the stochastic TDSLDA, with hopping 
between various time-dependent meanfields, and can ac- 
count for dissipation |53j . 

The most convenient quantity for studying the GDR 
is the strength distribution: S{E) = I (^^I-P'IO) - 
i?i^), where {v) are nuclear eigenstates corresponding to 
energies Ei,. In our case the operator is a sum of 
two operators depending on neutron and proton coor- 
dinates, respectively: Fr{r) = iV,- sin(k • rT-)/|k|, where 
T — n, p, |k| = 27r/L, L is a lattice size of the box, 
Nn = —Z, Np ^ N are neutron and proton numbers, re- 
spectively and A = N + Z. The operator Fr{r) generates 
all the odd multipoles, but the predominant contribution 
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comes from the dipole mode. In practice, depending on 
the symmetry of the ground state, we calculate the re- 
sponse of the nucleus to different external fields. Thus, 
for a spherical nucleus, we can chose Vp and r„ along any 
direction, for an axially deformed nucleus we compute 
two different responses with Vp and r„ along two sym- 
metry axes, while for a triaxial nucleus we compute the 
responses along three principal axes. Within our formal- 
ism, the external perturbation is added to the Hamilto- 
nian hT^„cr{r,t) => /iT-_o-CT(r, t) -|-FT-(r)/(i) and switched on 
adiabatically, where f{t) = C exp[—{t — 10)^/2] (time in 
units of fm/c) and C defines the intensity of the perturba- 
tion and which has to be kept sufficiently small to stay in 
the linear response regime, but large enough to excite the 
modes of interest. The amount of energy deposited into 
a nucleus was in the range 45-50 MeV. For such pertur- 
bations, which imply harmonicity of the excitation, the 
response function to the external perturbation F f(t) is 
given by j47]: S{uj) = Im{(5F(a;)/[7r/(a;)]}, where 5F[u) 
is the Fourier transform of the fluctuation of the expec- 
tation value: 5F[t) = {F)t - (F)o = / cPr5p{r,t)F{v). 
The extraction procedure of the strength function is sim- 
ple to state: obtain the self-consistent stationary solu- 
tion using the stationary SLDA solver to high precision 
so the spurious contributions due to global excitations 
(center-of-mass motion, rotations for deformed nuclei, 
proton and neutron pairing Goldstone modes) are de- 
coupled, and subsequently use this solution as input for 
the TDSLDA code where it is perturbed by the exter- 
nal field and evolved for a period of time, T. The en- 
ergy resolution due to taking the Fourier transform in 
the finite window T is 5E = 2T:h/T, and so it is im- 
portant to perform propagation to sufficiently large T 
(here T « 1600 fm/c and 5E « 0.8 MeV). We have used 
a cubic lattice box with L = 32.5 fm on a side and a 
lattice constant 1.25 fm, which allows us to obtain the 
ground-state energy with an accuracy of a few tens of 
kilo-electron-volts for a number of spherical light open- 
shell nuclei. The TDSLDA equations were integrated in 
time using a fifth order multistep method with a time 
step At « 0.12 fm/c chosen to maintain a 10~^ rela- 
tive accuracy, see Refs. ■ The changes in proton 
and neutron numbers during the time evolution change 
are \AZ\ < IQ-* and |AiV| < 10"^ respectively. TD- 
SLDA, when the external perturbation does not change 
the particle number, as in the GDR case, the particle 
number in the time-dependent solution is exactly con- 
served, and thus all the excited modes have exactly the 
same average proton and neutron particle number as the 
ground state. The method recently suggested in Ref. 
[17j . with the technical improvements described in Ref. 

, has the promise to become the leading approach for 
QRPA calculations in open-shell nuclei, because of the 
simplicity of its numerical implementation and its accu- 
racy. The solutions of the QRPA however are harder to 
interpret [2], since for a given nucleus these equations 



will provide solutions for AN = 0, AZ = 0, as well as 
pairing vibration type transitions with AA'^ = ±2 and/or 
AZ — ±2 respectively. TDSLDA has been applied so far 
to describe excitations with no change in particle number 
AA^ = due to time-dependent external fields which do 
not change the particle number [331 IM] , as well as excita- 
tions with AA^ ^ [3T] in the case of more complicated 
time-dependent external probes. 

From the TDSLDA solutions we construct the occupa- 
tion probabilities for both proton and neutron quasipar- 
ticle states: 




where k labels the proton and neutron quasiparticle wave 
functions respectively. Because of the complexity of the 
QRPA equations two types of approximations have been 
often used in the past in the literature. In one kind of 
approximation the proton and neutron densities are con- 
structed from 

Vk{r,<T,t) = ^Jnk{t)(f)t,{v,a,t), 

where 0fc(r, cr, t) satisfy TDHF-like equations. The occu- 
pation probabilities nk{t) are obtained either in a TD- 
BCS formalism, or are simply assumed to be frozen at 
their ground state values. One can easily show that in 
the case of TDBCS the continuity equation is not satis- 
fied [54j . and thus particle number conservation is vio- 
lated. The assumption that occupation probabilities are 
frozen to their ground-state values is violated strongly as 
well, as the two plots in Fig. [ijhere amply demonstrate. 
A typical example of this frozen occupation probabilities 
approximate approach is in Ref. .55], in which one can 
find a calculation of the GDR in the triaxial nucleus ^®*0s 
with SLy6, which is very similar to SLy4. The shape of 
the GDR line shape calculated in Ref. [55] is noticeably 
different from both the experimental one and our results; 
see Fig. [2] We notice that there are experimental indica- 
tions that a pygmy resonance exists in the nucleus ^^^Yb 
[56) . We did not find any evidence of a pygmy resonance 
using either of the two Skyrme energy-density functional 
parametrizations SLy4 and SkP, and neither did Terasaki 
and Engel using SkM* [T^]. In Ref. [T^] pairing corre- 
lations were treated using a sharp cutoff, a prescription 
that shows a strong dependence on the cutoff energy, at 
least in the case of low-lying states [57]. The approach 
adopted by us is free of such cutoff dependency in the 
pairing channel, once the cutoff energy is chosen roughly 
above 50 MeV pi". 

In Fig. [2] we show the results of our calculations and 
compare them with experimental data [5S]. For the axi- 
ally deformed nuclei, we performed two calculations per- 
turbing the systems along the longer and the shorter axis, 
respectively. The triaxial nucleus -"^^^Os required three 
runs accordingly. The experimental data include sev- 
eral effects: the coupling of the collective strength to the 
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FIG. 1. (Color online) The time-dependent proton and neu- 
tron occupation probabilities of a set of quasiparticle states 
around the Fermi level for ^•^^^U calculated as described in 
the main text with an approximate energy density functional 
using the SLy4 parametrization. 



Ip — Ih states, especially to the low-lying collective vi- 
brations, Landau damping (F^ « Hvp/R, where vp is 
the Fermi velocity and R the nuclear radius), particle 
escape into the continuum (F^), and coupling to more 
complex states (F^). The theoretical Thomas-Reiche- 
Kuhn sum rule for these nuclei are: 2.49 (^^^Yb), 2.81 
(is^Os), and 3.38 (^^^U) in MeV-barn respectively [5]. 
The corresponding energy integrated photo-absorption 
cross-sections in the interval 8- to 20-MeV are: 2.47, 2.73, 
and 3.06 in experiment, 1.94, 2.26, and 2.49 for SkP, and 
1.84, 2.01, and 2.35 for SLy4 respectively. Both SkP and 
SLy4 forces underestimate the energy integrated cross- 
section in this energy interval, even though both Skyrme 
parametrizations adequately reproduce the average peak 
position and even the width, in spite of great difference 
in the isoscalar nucleon effective mass m*^ = m and 
m*^ = 0.7 m, but very similar isovector effective mass. 




12 16 20 
E [MeV] 



12 16 20 

E [MeV] 



TO,-, 



0.74 m and m*^ = 0.8 m for SKP and SLy4 respec- 



FIG. 2. (Color online) The calculated photo-absorption 
cross-section (solid black line), using two Skyrme force 
parametrizations for three deformed open-shell nuclei and the 
experimental (7, n) cross-sections (solid purple circles with 
error bars), extracted from Ref. [58]. With dashed (green), 
dotted (red) and dotted-dashed (blue) lines we display the 
contribution to the cross-section arising from exciting the cor- 
responding nucleus along the long axis, the short axis (mul- 
tiplied by 2 for the prolate nuclei ^^^Yb and '^^^U) and the 
third middle axis in the case of the triaxial nucleus ^**0s. 



lively. The use of periodic boundary conditions forced us 
to use |k| = 2ti/L in Fr{x) = A^t- sin(k • r7-)/|k|, instead 
of |k| — > 0, which leads to an underestimate of the cross- 
sections by w 25%. With this trivial correction (obtained 
by calculating the energy weighted sum rule for the ex- 
ternal probe we obtain 2.63(6) (^^^Yb), 3.11(7) {^^^Os), 
and 3.55(8) (^^^U) for SkP and 2.50(5) (^^^Yb), 2.77(5) 
(is^Os), and 3.35(8) (^^Sy), with the SLy4 results being 
in slightly better agreement with experimental results. 
The mixing with more complex configurations (F"^), will 
lead to some depletion of the transition strength out of 
this energy interval. 
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